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I. INTRODUCTION 



Nano plasmas can now be readily produced in laser irradiated clusters, and new physical phenomena have come 
into focus experimentally as well as theoretically. Interactions between laser fields of 10^'^ — 10^^ W cm~^ and clusters 
have been investigated over the last few years, see Refs. [H-@. After laser interaction, extremely large absorption 
rates of near ly 100%, see as well as x-ray radiation, see Refs. [ii|-[i3j were found. In pump-probe experiments, 
e.g. by Doppner et al. ^ and Fennel et al. [1^, the absorption rate of a second laser pulse is strongly dependent on 
the time delay what is caused by the dynamical properties of the expanding cluster. We will discuss the dynamical 
response function of the electrons in a nano plasma that is responsible for scattering and absorption of electromagnetic 
radiation. 

Collective electronic excitations of the nano plasma usually are interpreted as Mie resonances of a homogeneously 
charged sphere. Absorption cross section experiments by Xia et al. [l9| show multiple resonance structures indeed. 
The effect of collective electron motion can also be seen in ultraviolet (UPS) and x-ray photoelectron spectroscopy 
(XPS) experiments, see [lO], which are used to detect binding energies of core level electrons in small metal clusters 
|2ll [22j . In fusion related experiments by Ditmire et al. [IJ , Grillon et al. [13] , as well as Madison et al. [1^ ignition 
processes are started via irradiation of deuterium clusters. Collective effects in the optical response are discussed in 
the context of metallic nanoshells by Hoflich et al. [13 as well as nanocavities by Maier et al. [27|. 

In theoretical calculations of finite systems, see Raitza et al. l28l-|3ll. a more complex resonance structure was 
found. Earlier investigations by Reinhard et al. [s^l and Kull et al. |33l. fs^] led to comparable results. The method 
of resonance structure analysis using spherical harmonics is known from the discussion of g iant dipole resonances of 



nuclei, see Reinhard et al. [35]. Quantum and semi-classical methods, see Refs. [36|-[38[, respectively, were used 
to investigate the cluster excitation via laser fields. Collisional absorption processes in nano plasmas have been the 
subject of theoretical investigations by Hilse et al. [s^. Usin g d ensity functional theory (DFT) calculations, the 
electronic structures of cold clusters were analyzed by Ekardt |40|. Kiimmel et al. HJ], Brack et al. as well 

as Krotscheck et al. [i^. The damping of collective electron oscillations was investigated by Ramunno et al. (i^ ] 
emphasizing the importance of collisional processes beside the Landau damping. 

In this work, molecular dynamics (MD) simulations will be used to study nano plasmas in metal clusters. Clusters 
consisting of 55 up to 1000 sodium like atoms are considered after short pulse laser irradiation with intensities in 
the order of 10^^ Wcm^^. Properties of the nano plasma are mainly determined by the dynamics of electrons which 
are bound to the cluster but ionized from the former atoms, comparable to conduction electrons in bulk systems. 
As already shown in earlier publications, see f28j , plasma parameters as known from bulk (temperature and particle 
density) but also the cluster size and net charge are justified for characterization since the electrons can be assumed to 
be in local thermal equilibrium within time scales considered here. We focus on parameter ranges where the plasma 
can be treated classically. Strong correlations are taken into account via collisions of all particles. Concepts that have 
been well established for infinite bulk systems near thermodynamic equilibrium have to be modified for applications 
to finite systems, e.g. clusters. In particular, we are interested in the dynamical structure factor and the response 
function for such finite nano plasmas. In order to bridge from finite systems to bulk plasmas, we investigate size effects, 
e.g. in the dynamical collision frequency. First results in this direction have been reported in Refs. (28l. [isl. l46j . 

In Sec. II, correlation functions and their relation to optical properties of homogeneous bulk plasma are introduced 
as far as it will be of interest to extend the approach to finite systems. Expressions will also be used for comparison 
with nano plasmas in the limit of large clusters. Sec. Ill explains the restricted molecular dynamics (RMD) scheme 
for the calculation of the particle trajectories from which the total and bi-local current density correlation functions 
are determined. Symmetries in the correlation matrix discussed in Sec. IV can be used for an improved statistics. 
The decomposition of the correlation matrix into eigenvectors and eigenvalues is interpreted as a decomposition into 
collective excitation modes. In Sec. V, the excitation modes will be characterized with respect to spherical harmonics. 
In the further analysis, we focus on modes with a dipole moment, which are also seen in the total current density 
auto-correlation function. First results for resonance frequencies and damping are presented. Regarding the dipole- 
like modes, the spatial structure at the selected resonance frequency will be discussed in Subsec. A and Subsec. B. 
Conclusion and outlook are given in Sec. VI. 



II. LINEAR RESPONSE THEORY OF PLASMAS IN EQUILIBRIUM 

Within linear response theory as derived by Kubo et al., see [U I3, the reaction of a many-particle system 
to weak external perturbations can be related to the dynamical behavior of fluctuations in thermal equilibrium. 
Denoting the equilibrium statistical operator with po , we introduce the two-time correlation function of the fluctuations 
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6Ai{t) = Ai{t) — Ti{AiPQ} as the Kubo scalar product 

{A,{t)-AM)^P I A\Tv{5A,{t)5Aj{\h(3\)po), (1) 
Jo 

where the time dependence is given in the Heisenberg picture. The indices i and j identify quantum observables. 
In particular we consider local properties so that they contain also the position r. In case of i = j, it is called 
auto-correlation function (ACF). 

In the classical case, equilibrium two-time correlation functions can be calculated according to 



1 r 

(A,(f,t);A,(r',0))= hm - / dr M,(r, r + i) • (r', r), 



(2) 



where we assumed ergodic systems - the ensemble average can be replaced by a time average. The spectrum of the 
equilibrium correlation function (Aj(r); Aj(f*))^ then results from Laplace transformation. 

We consider an induced electron density fluctuation Snc{f, t) = ne{f, t) — nefi{r) at time t as the deviation from the 
equilibrium density distribution nefi{f^ due to an external potential J7ext(r*,t') at times t' < t. Close to equilibrium, 
the correlation between the external potential and the induced density fluctuation is only dependent on the time 
difference Ai = t — t'. Thus, one is able to discuss its spectrum after Laplace transform. In the same way, the induced 
electrical current density jo{f,t) is related to the external electric field E{r' ^t'). Via Kubo's theory, these induced 
quantities, S{nc{r))^ and S{j{r))i^, can be expressed within linear response, see [ll], as 



S{n,ir))^ = /3y dV (<5ne(f);fee(r'))^ C/ext(r', w), (3) 
= /5 / (J(r); jV))^ E{t^,u:). (4) 



The spectrum of the density fluctuation correlation {Snc{r); Shc{f^))uj is related to a scalar response function. The 
current-density correlation {j{r)] j{f'))^ represents in general a tensor due to the directions of the current density 
vector. 

Before considering non-local response functions, we shortly mention homogeneous systems. Properties of the bulk 
plasmas with electron density ric and inverse temperature /3 are only dependent on the difference of the positions 
Af = r ^ f" . Thus, after Fourier transform of the spatial difference Ar, the correlations are dependent on a wave 
vector k. 

The dynamical structure factor is directly related to the density fluctuation correlation, as 

S{k,L^) ^ ^^j^{5nk\5nk)^ (5) 

with N the number of particles. For further relations to the dielectric function and the optical response of a homo- 
geneous bulk plasma see [i^. Note that the density fluctuations Eq. ([3]) as well as the density correlation function 
in Eq. ([5]) can be expressed in terms of the current-density correlation function via partial integration and using the 
continuity equation. Thus, the dynamical structure factor is divided into a static part So{k) and a dynamical part 
which is directly related to the longitudinal part of the current-density correlation function 

It is of fundamental interest to describe the collective behavior of the system as response to external fields, in 
particular emission, absorption and scattering of light. In bulk systems, the wave vector and frequency dependent 
response function reads 

fc^ -II -I 

which can be evaluated using quantum statistical approaches such as Green function theory, see [i^, or numerical 
approaches such as MD simulations, see [13] ■ As collisions are relevant in stron gly correlated systems, the dynamical 
collision frequency ^{ijj) is derived and appears in a generalized Drude formula |5lL Is^ 

lim x(fc, c.) - . (8) 



^a;2 - Lo^^ + iuju{uj) 
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In the classical case, the current-density correlation function has been extensively discussed in the long wavelength 
limit /c — >■ applying MD simulations and perturbative approaches. Exemplarily, we refer to (soj . 

The state of a homogeneous one-component plasma in thermodynamic equilibrium is characterized by the nonide- 
ality parameter F = e^(47rnc/3)^/'^(47reo^B7c)~"'^ and the degeneracy parameter Q = 2mekBTeh~^ {3Tr^nc)~^^^ , Tc is 
the temperature of the electrons. Considering the response function x{k, oj) in the long wavelength limit, a sharp peak 
arises at the plasmon frequency Wpi, see Eq. ([5]). For finite wavelengths, the resonance is shifted and can be approxi- 
mated by the so called Gross-Bohm plasmon dispersion for small wave numbers k, see [H,[53|, uiik) « uj-p\ + 'ik'^ / + ... 
with the Debye screening length = [nee^/feofcaTe)]""'"^^. This relation has recently been revisited with respect to 
the relevance of collisions by Thiele et al. (S^. According to Eq. ([5]), the general behavior of the response function 
x{k,uj) in the long-wavelength limit is closely related to the collision frequency which is relevant in non-ideal plasmas, 
see [53 . [ssj . In the two-component plasma, a phonon mode can arise in addition to the plasmon excitations [scj . 

The response function x(fc, i^) and the related dynamical structure factor S(k,uj) as well as the optical properties 
have been intensively investigated for electron-ion bulk systems, see Refs. [ssl ISTj. In this work, the inhomogeneous 
case of finite clusters in local thermal equilibrium will be discussed. The response of inhomogeneous systems is not 
only dependent on the difference of the positions, but on f and f* separately. Therefore, spatially resolved current 
density correlation functions {j{r); j{r'))^ can not be diagonalized by spatial Fourier transform. Instead of plane 
waves, other basis functions have to be found in order to characterize the collective excitations of electrons. 



III. MD SIMULATIONS OF FINITE PLASMAS 



Finite plasma systems have been investigated using the restricted molecular dynamics (RMD) simulations, see 
Raitza et al. [28]. A two-component system of singly charged ions and electrons will be described using an error 
fimction pseudo potential for the interaction of particles i and j 

where Zi is the charge of the ith particle. The Coulomb interaction is modified at short distances, assuming a Gaussian 
wave function for electrons motivated by the account of quantum effects. Considering a sodium like system, the 
potential parameter A — 0.318 nm was chosen in order to reproduce the ionization energy of Ip — Vci{r ^ 0) = —5.1 eV 
for solid sodium, as already discussed for MD simulations by Suraud et al. [ssj . 

The velocity Verlet algorithm [s^ was applied to solve the classical equations of motion for electrons and ions. This 
method takes into account the conservation of the total energy of the finite system, as long as there is no external 
potential. To follow the fast electron dynamics, time steps of 0.01 fs were taken to calculate the time evolution. 
Contrary to bulk MD simulations no periodic boundary conditions are applied. 

Icosahedral arrangements of 55, 147, and 309 ions, see [l^, were considered as initial configuration for the ion 
positions. For these nearly spherically, homogeneously distributed ions, the ion density typical for solid sodium is 
given by an ionic next neighbor distance of do = 0.212 nm. In addition, randomly distributed ion configurations 
within a given sphere were considered for comparison and the number of ions was increased up to 1000 particles. 
Starting with a neutral cluster, the electrons have been positioned nearby the ions with small, randomly distributed 
deviations from the ion positions. 

To simulate experiments where clusters are excited by short pulse lasers, MD simulations are performed under 
the infiuence of an electric field, assuming a Gaussian shape and pulse duration of about 100 fs. Due to the largely 
increased kinetic energy of the electrons, ionization processes occur. After the laser field is switched off, the ionization 
degree of the cluster is determined by the number of electrons found outside the cluster radius with positive total 
energy, so that they can escape from the cluster. Due to ion excitation on larger time scales, a slow expansion of the 
positively charged cluster is obseved , leading to Coulomb explosion experimentally. 

Considering the single-time properties, it was found in ^28^ that already local thermodynamic equilibrium (LTE) 
is established within a few fs after the electron heating. In particular, at each time step, the momentum distribution 
of electrons is well described by a Maxwell distribution, and the spatial density profile agrees with a Boltzmann 
distribution with respect to the average potential that is determined by the actual ion configuration and the self- 
consistent electronic mean field. The fact that electrons are considered within sub fs time intervals, while the ion 
configuration remains nearly unchanged, enables us to separate the electron dynamics from the ion dynamics. 

Subsequently, the dynamical properties of the electron subsystem can be calculated for a frozen ionic configuration 
thus referring to a specific time. This is considered as an adiabatic approximation to the true dynamical properties 
of the electron subsystem which have to take into account the slow change in the ion configuration. More rigorously, 
non-stationary time dependent correlation functions have to be treated for the full charged particle system. 
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FIG. 1. (Color online) Simulation data of the cluster charge Z (symbols) and a power fit (solid line) for clusters with different 
numbers of ions M at different temperatures is shown. 



Using the RMD simulations scheme as introduced in [28|, the ions are kept fixed acting as external trap potential. 
Starting from an initial state, the many-electron trajectory {n(t),pi(t)} is calculated, solving the classical equations 
of motion of the electrons. From this, all further physical properties of the electron subsystem inside the cluster 
are determined. Within RMD simulations, we consider no temporal variation of the plasma parameters that are 
determined by the frozen ion distribution, the electron temperature and the degree of ionization. A long-time run can 
be performed in order to replace the ensemble average by a temporal average. This has been successfully done for the 
single-time properties such as the momentum distribution and the density profile, see (281] and will now be applied to 
the two-time correlation functions. 

Using classical MD simulation techniques, the results are valid for non-degenerate plasmas. This restricts the 
temperature range to T > 1 eV where our simulations can be compared with realistic sodium clusters. Values for the 
plasma parameter F > 1 can be treated since we are not confined to the weak coupling limit as, e.g., in perturbation 
theory. 

In our RMD calculations, we start from a homogeneous ion configuration (icosahedral or randomly distributed) 
inside the cluster at fixed ion density. In the case of random distribution, we perform averaging over different initial 
configurations of ions. The Langevin thermostat was used to heat the electrons at an initial stage. We have chosen the 
Langevin thermostat introducing a friction term with suitable sign to adjust the intended kinetic energy. Furthermore, 
a random source term is applied that thermalizes the system. Hot electrons are emitted during this stage so that 
the cluster becomes ionized. Evaluating the trajectories of electrons, sufficient time of about 200 fs has to be allowed 
before a stationary ionization degree is established. Then the thermostat is switched off and a data taking is performed 
using an ensemble at fixed number of density, volume and energy. It is checked that the mean cluster charge Z and 
the system temperature do not change any more. In Fig. [TJ the cluster charge Z depending on cluster size A'i is 
shown. A power fit Z(iVi) = AN^ with, for example, A = 0.165 and B = 0.197 for Te =1 eV shows the trend of the 
size dependent ionization degree. 

Using the trajectories of all electrons obtained from the RMD simulations scheme, the local current density 
je {fi t) at position r was calculated for each time step t 



which is the sum over all electron momenta pi inside a small volume AVp at position r, where <5AVp(n(i)) = 1, and 
SAVri'Tiit)) ~ for electrons found outside AVf. The size of the volume determines the spatial resolution of the local 
current density j(.{r,t). However, it must be taken sufficiently large to reduce statistical fluctuations. 




(10) 
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The bi-local correlation tensor of the normahzed spatially-resolved current density is calculated according to Eq. ([2]) 



as 



(je(r,0);je(r',t) 



J2^=i io{r, i ■ t) j'^{r',i ■ T + t) 



Nr ■ if.) 

with jc the total current density. Typical values are — 10^ — 10^ and r ^ 0.1 is. Its Laplace transform reads 



(11) 



(Jc(r);Jc(r')).= 



dti 



In the following, we restrict ourselves to the diagonal components (jc ; Jc)lj of this tensor, where only parallel compo- 
nents of the current density vectors are correlated as already introduced in Sec. II. As it will be shown in the following 
sections, this bi-local current-density correlation is important to understand the excitation modes of nano plasmas. 
The non-diagonal components of the bi-local correlation tensor are small in comparison to the diagonal components. 
Beside the bi-local current density correlation function considered here, the bi-local density fluctuation correlation 
{5n{r),5n{f'))uj as well as the bi-local force correlation {F{f),F{r'))^ are useful quantities in the context of optical 
properties. These correlations can be evaluated from the trajectory in a similar way and are related to the bi-local 
current density correlation. This will be discussed in an upcoming paper. 

Because of the spherical symmetry of the cluster geometry during excitation and expansion, the volume is divided 
into sections AVf? according to N^., Ng, equidistant intervals of spherical coordinates, i.e. the distance r to the 
center of the cluster, the inclination angle as well as the azimuthal angle </>, respectively. The cluster radius Ri is 
given by the root mean square radius of ions according to = 5/3(r^). The sections are numbered by a single counter 
a = NfpNe (fc — 1) + {j — 1) + i with three independent counters according to the three coordinates: i = l.-Ncj,, 

j — l..Ng and k — l..Nr. With respect to Eq. (|lip the bi-local correlation matrix Da-a'{t) = {jc{f'a,^)]jc{f'a',t)^ for 

the spatially resolved cluster and its Laplace transform Da;a'{oj) = /q°° dte"^* Da-a'{t) have been calculated. 

The total current density ACF can be calculated from the trajectories directly. Please note, that it can be also 
calculated from the bi-local current density correlation matrix 



(12) 



,j',k' 



(13) 



using the cluster volume V^i = ^-Rf and the individual cell volumes 



27r(-i-l)/Af^ 

27r f Ri 



TTj/Ne 



Rik/Nr- 



AO I dr sin 9 



3iV0 \Nr 



(3fc2 - 3fc + 1) cos 



Ne 



{3-1) 



(14) 



The consistency of these expressions has been checked throughout our explicit calculations. 



IV. FROM BI-LOCAL CORRELATION FUNCTION TO EXCITATION MODES 

In the following, we discuss calculations for the current-density ACF, Eq. ([T^. and the bi-local current-density 
correlation spectrum Da^a'{^)- Exemplarily, we present results for the Naaog cluster at electron temperature — 1 
eV, cluster charge Z — 16 and ionic density ni — 2.80 • 10^^ cm~'^. The electrons form a nano plasma with nonideality 
parameter F — 6.964 and degeneracy parameter O = 0.664. Starting with a solid density cluster, these are typical 
parameters obtained directly after the interaction with a short pulse laser of 100 fs duration and intensity of / = 5-10^^ 
Wcm~^. Calculations of other cluster sizes will be presented in the following sections. 

The real part of the total current-density ACF Re(ji';ji')^ is shown in Fig. [21 Three maxima are obtained. This 
feature differs from the bulk behavior and is interpreted as different resonances of the electron system. To investigate 
the origin of the different maxima as collective excitations of the nano plasma, the bi-local current-density correlation 
matrix was calculated as well. 

The following spatial symmetries in the matrix Da:a'{^) were found 

Di^j^k; i' ,k'i^) — + |i-i'| + lj-',fc'('^)j (15) 

Di^j^k- i'J',k'{^) — Di,Ne-j + l,k- i' ,Ng-j' + l,k' {i^) , (16) 
Dij^k; i'j',fe'(w) = Di,j,k'; i'j',fc(w). (17) 
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FIG. 2. (Color online) Frequency spectrum of the real part of the current-density ACF of a Nasog cluster at electron temperature 
Te = 1 eV, cluster charge Z = 16 and ionic density ni = 2.80 • 10^^ cm~^. 

In our case, the N^^^ elements of the full matrix can be reduced to A'i„d — {Nr Ng + 1) Nr Ng {N^ + Af0mod2) /4 
independent elements due to the symmetries Eq. (jlSp - Eq. (I17|) . thus improving statistics via averaging equal elements. 

Because of the different size of section volumes in spherical coordinates there are large variations in the mean 
number of particles in a section. Provided that we have A^e = 50 — 1000 electrons and A'scc = Ncf,NeNr = 128 sections 
the average number of particles in some sections can be even smaller than unity. In this case, the local current density 
Eq. ([TU|) is affected by strong fluctuations due to the discrete number of particles. This problem is reduced, when you 
consider the current Jc{r, t) — jc{r, t)AVp as the contribution of smaller cells will be damped. Therefore, we used the 
non-normalized form of the correlation function for further analysis 

for which the following eigenproblem was solved, 

Reifa,a'(w)*M,a'(^) = K f,{uj)^ f,^a{uj) ■ (19) 

a' 

Thus, the matrix is decomposed into eigenvectors 5'p.a(aj) = i_i,ij.k{^) as well as their eigenvalues K^{uj) at each 
frequency. The eigenvectors represent the spatial structure of the mode (^;^,i.j,fc(w) — > ^^(r, w)). The orthonormality 
condition 

j dV-*^ (r , Lo) (r , 0.) = 5^,^' . (20) 

holds. 

For two selected frequencies, the 10 strongest eigenvalues of the Naaog cluster are shown in Fig. [31 At cj = 4.42 
fs""'^ (black), a resonance frequency was found with one outstanding, leading eigenvalue. The second and third largest 
eigenvalue are of same strength, which suggests degeneracy due to the symmetry of the correlation matrix. At off- 
resonant frequencies, i.e. at w = 5.50 fs^^ (shaded, red online), all eigenvalues are of the same order of magnitude. 

In Fig.|4](a) the strongest eigenvalues K^{uj) of the Naaog cluster are shown in dependence of frequency. They are 
colored according to their strength and numbered ascending with descending strength. In the shown frequency range, 
modes Kfj_{uj) with well defined maxima are found. The spatial oscillation structure can be identified by analyzing 
the eigenvectors. 

In Fig. |4] (b) , the spectra of eigenvalues are sorted in an alternative way, according to the spatial structure of the 
eigenvector which is obtained over the whole frequency range. Overall, the black solid mode is the strongest. Its 
resonance frequencies are also found in the total current-density ACF (indicated via vertical blue dashed lines) and 
are therefore of particular interest. Resonances in the total current-density ACF, shown in Fig. [51 are only possible 
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FIG. 3. (Color online) Eigenvalues of the Naaog cluster sorted by size for the resonant case at lj = 4.42 fs~^ (black) and a 
non-resonant case at cj = 5.50 fs~^ (shaded, red online), for parameters as in Fig. [2] 



■ CO = 4.42 fs 
• co = 5.5fs"^ 




in the case of non-zero total current, which is caused by a dipole-like oscillation. Thus, resonances which are seen in 
the total current-density ACF are oscillation modes with a dipole moment. Other resonance structures, for example, 
are breathing modes that have no dipole moment. After characterization of the resonance structures, the dipole-like 
resonances will be investigated in more detail. 



0.004 I ^ ^ , ^ , ^ 




CO [fs"'] 0) [fs"'] 



FIG. 4. (Color online) Spectrum of the 6 highest eigenvalues A'^(tj) of the Nasog cluster for same parameters as in Fig. [2] (a). 
Eigenvalues of the same cluster selected in terms of the corresponding spherical harmonics Yi,m(^, 0) (b). 



V. ANALYSIS OF THE COLLECTIVE MODES 



The decomposition of the locally resolved current correlation matrix into eigenvalues i4r^(a;), as shown in Fig. UJ 
gives a very complex set of resonance structures in comparison to the ID case, see [13] . The spatial mode structures in 
ID chains were characterized by their wave number k. To analyze the more complicated spatial oscillation structure 
of 3D clusters, a spherical Fourier decomposition of the eigenvectors into the spherical Bessel function ji{kn^ir) and 
spherical harmonics Yi^„i{(^, (f) performed according to 



JV„ Ni I 
Ti=l i=0 m= — l 



(21) 
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where S'n,!,m(w) is the spherical Fourier component with ordinal numbers n^l^m. The normalization factor Nn^i as 
well as the wave number fc„ ^ are chosen in the way that the eigenvector has a root at the cluster surface. 

In Fig. |4] (b), the four strongest eigenvalue modes are characterized by pairs of ordinal numbers I, m which deter- 
mine the main angular part of the eigenvector by the spherical harmonics yi^m(0, 0). The leading dipole-like mode, 
represented via solid black lines in Fig. |4] (b), is characterized by the overlap of the spherical harmonic functions 
Yo,o{(^i4') and Y2fi{0,(p). For the Naaog cluster, one can find three resonance frequencies which are identical to the 
ones found in the total current-density ACF. The latter are indicated by vertical dashed lines (blue online) in Fig. |3] 
(b). 

In our investigations, we looked at other cluster parameters as well and found similar behavior. Comparisons will 
be made in the following chapters. For further analysis of the exication modes, we now consider a larger cluster 
consisting of 1000 ions. There, four pronounced dipole-like resonances were found. In Fig. [Sj the spatial structures of 
the current-density jc {r) ^ AV(r) '^^ shown for the Naigoo cluster at the resonance frequencies of the leading dipole-like 
mode. The behavior is shown in the z — a;— plane at a fixed azimuthal angle (f) on which it does not depend. 



w = 4.80 fs~^ oj = 6.18 fs~^ uj = 9.15 fs~^ lo = 8.23 fs~^ 
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FIG. 5. (Color online) Selected eigenvectors of the dipol-like mode in the Naiooo clusters for same parameters as in Fig. [2] but 
Z = 19. 

At the resonance frequency = 4.80 fs, the electrons are oscillating with a current density jc{r) = v(r)nc{r). As 
shown in Fig. [5] (a), all electrons of this mode are moving in the same direction and no nodes can be seen. Assuming 
a constant velocity field amplitude v = const, the change of the current density with distance r is directly related to 
the density profile riair) of the electrons. 

The modes in Fig. [5] (c) and (d) are similar to a plane wave oscillation of electrons, but trapped inside the cluster. 
To identify a wave number of the plane wave oscillation, a Fourier decomposition of plane waves in ^-direction was 
done. A maximum at fc = 1.6 nm~^ and k — 4.7 nm~^, respectively, is found which identify the wavelengths of the 
plane wave oscillations. Only in the large cluster with 1000 ions, a plane wave oscillation with higher wavenumber 
was found. All other modes can be seen in smaller clusters as well. The resonance structure in Fig. [5] (b) looks like a 
mix of the first and the third resonance structure. 

We want to point out one further feature of the mode spectra in Fig. |4] (b). The dashed red line represents in 
fact two resonance structures with exactly the same eigenvalues at all frequencies. The eigenvectors are orthogonal 
since they are characterized by the same spherical harmonic function Yi^i(0, 0) but have a phase shift in 0-direction: 

(6», 0) = yI"^} (6*, <^ f ). Further degenerations are obtained for weaker eigenvalue modes as well. 
All eigenvectors \E'^(r, w) are decomposed into a superposition of spherical Bessel functions ji{kn^ir) with a set of 
ordinal numbers n. No leading ordinal number n was found, which characterizes the spatial resonance structure in r 
direction. 



A. Resonance frequency of the rigid oscillation 

The total current density ACF shown in Fig. [5] as well as the leading eigenvalue mode in Fig. 2] (right) show 
the strongest resonance at the frequency cjr k, 4.42 fs^^. This resonance belongs to the dipole-like mode with the 
eigenvector shown in Fig. [S] on the left hand side. We will now analyze this collective excitation mode in terms of a 
rigid oscillation. 

The electrons with density profile nc(r) are assumed to move nearly rigidly in the external potential VcyLt^ciir) due 
to the fixed ions. The potential energy of the electrons due to a small shift with respect to the ions reads 

t/e(z) = / dV7ie(f) K=xt,ci(r - ze,). (22) 
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The change of the potential energy Uc{z) in z direction is due to the restoring force on the electron profile. In harmonic 
approximation of the equation of motion, the resonance frequency is identified as 



52 U,{z) 



dz^ 



(23) 



z=0 



For small rigid shifts z — > 0, assuming radially dependent electron density profiles and external potentials in Eq. 
the integration over the angular dependence of the potential energy calculation can be executed. The resonance 
frequency Eq. (j23p is then given according to 



47r 



(24) 



As a first example for a density profile, we assume a homogeneously charged ion sphere with radius i?i = 
(3Afi/(47rni))^/'^ and an electron sphere with radius Rc — (3Afo/(47rnc))^^'^- The densities of the electron and 
ion spheres is equal (rio = ni). Therefore, the difference of ion and electron radius is determined by the cluster charge, 
basically the difference of the simulated electron number and ion number Ni. Thus, in the case of positively charged 
clusters, as discussed here, the electron sphere radius is smaller than the ion radius {Rc < Ri). The error function 
potential Eq. Q was taken as electron-ion-interaction potential for the calculation of the resonance frequency, as it 
was used for the MD simulation as well. The resonance frequency than reads 



uj^{Ri,Rc 



Rf + i?e r -^i + 

■erf 



2i?3 



Rf — Rc r f R'l ^ Rc 
-ert 



2i?3 



Rn 



sinh 



2i?i_R(. 



Ai?ii?cCosh 



/ 2R;R„ 

— 



(25) 



In the limit of large clusters with high number of ions the resonance frequency equals the Mie frequency, 
limfl;_j.tx, ajR (i?i) = WMio- Assuming only a weak charged cluster, the sphere radii have nearly the same size (i?o — > ^^i) 
and the system is nearly neutral. The limit for small clusters, down to just one atom, depends strongly on the pseu- 
dopotential. In our case, the resonance frequency limjvj^i i^R(A'i) = e'^/Aneo ■ 4/ {3^/^^X^mc) is due to the oscillation 
of a single electron in the ionic error- function pseudo-potential Eq. 

In Fig. ini (a), the resonance frequency wr of the dipole-like mode is shown in dependence on the size of the ion 
sphere. Results from MD simulations (empty circles) for Nass, Naaog and Naiooo cluster at Ui 



Eq. for ion densities of 



-)22 

2.15-1022 



well as the Na55 cluster at rii — 2.15 • 10 cm are shown. The resonance frequencies have been calculated using 



cm 3 (solid shaded line, red online) and ni = 2.80 • 10 cm ^ (solid black 



line). The limits of large clusters, the Mie frequency uj 
to the two densities. 



2 
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e ni/{3eomc), are given as dotted lines colored according 
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FIG. 6. (Color online) Cluster size dependent resonance frequency WR,(_Ri) (a). Simulation results (empty circles) and analytical 



calculations (solid lines) using Eq. (|25|) are shown for ni 
(black). The Mie frequencies are given as dashed lines. 



2.15 ■ 10^^ cm"^ (shaded, red online) and m = 2.80 • 10^^ cm- 
Numerical calculations using Eq. (|24p are presented (full dots) 



Dispersion Eq. ((32]) of a plane wave in a homogeneously charged sphere (solid line) for — 2.1 
simulation results (empty symbols) and bulk plasmon frequency (b). 



10 cm as well as 



Additionally, the electron density profile ^^(r) was deducted from MD simulations for all cluster sizes and used to 
derive the resonance frequency solving Eq. (j24[) numerically. As a result (full circles in Fig. [5] (a)), the resonance 
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frequency of the dipole-like mode is obtained with a deviation to the direct simulation results of less than 5%. Using 
homogeneously charged ion and electron spheres leads to reasonable agreement in the limits of large clusters as well 
as for small clusters. Taking the spatial structure of the density profile into account, there is good agreement with 
the direct simulation results in the intermediate cluster size regime as well. 



B. Dispersion of the plane wave mode 

While the Mie-like resonance, discussed in the previous subsection, is almost spherically symmetric, we obtain an 
increasing plane wave character of the dipole-like mode with increasing frequency. The third resonance frequency of 
the total current density ACF for the Naaog cluster at wr, = 8.14 fs~^, see Fig.[2l is mainly caused by a plane wave like 
eigenvector, which is similar to the eigenvector of the Naiooo cluster, shown in Fig. [5] (right). This mode is discussed 
by Kresin et al. as compressional volume plasmon. Here, oscillations of electrons in opposite directions must 
be taken into account for the analytical calculation of the resonance frequency. We assume homogeneously charged 
spheres for the electrons with radius i?c and for the ions with radius R\ as it was already discussed in the previous 
subsection. The electron motion is treated as a hydrodynamical liquid using the Euler equation 

djir.t) ,. r^,^ , ^.^ ,1 1 , X nJfA) 
' ' = -div "-^^ +^ ^ +^ „,.„j„/^+A ' ' 



dt 



i{f,t)®v{f,t) gradp(f,t) — —gT:a,AVc-Kt{r,t), (26) 



where j{f,t) — nc{r^t)v(f,t) is the spatially resolved current density of the electrons, p{r,t) is the pressure of the 
electron gas and Vext — K:xt,ci + K^xt.cc is the external potential, composed of contributions from the electrons and 
ions. Using the following ansatz 

j,(r,<) = <5jeV("), (27) 

v,if,t) ^ 6ve,&^'''-'^'\ (28) 

(r , t ) = ne,o (r) + SnJ<^'''~'''^ , (29) 

Kxt (r, t) = Uext,o + SV,^t (r), (30) 

we consider small perturbations in z-direction restricting ourselves to longitudinal effects. One is able to linearize the 
Euler equation. The system is assumed to be in LTE, described by the quantities nc,o(?'), jo(^ = 0, vo{r) = as 
well as Vext,o('")- Electrons are moving in the external field of ions and in the mean field of electrons. The external 
potential is 

/ ^ "c,o(ri)\ . ,^ ,1 /".g^ 



14xt(r,t) = J d'n l^n,{n) - ^"y^J K=,i(fi - ^) + 3 _/ d^ridn,{n,t)V,Ari - r), 

K=xt(r, t) = yext.o(0 + 6V,^t{r, t). (31) 

The external potential has a equilibrium part and a perturbative part 514x1(^,0: which is mainly dependent on the 
linear density perturbation 5ni^{f,t). 

Assuming Boltzmann distribution we express the ideal gas pressure of the electrons p{f, t) via the electron density. 
Using the equation of continuity, one is able to express the Euler equation in terms of linear perturbations of the 
density. Thus, the equilibrium part of the external potential compensates the pressure term on the right hand side of 
Eq. (|26|) . Restricting ourselves to linear perturbations of the Euler equation, only the third term on the right hand 
side of Eq. (pS)) remains, which is connected to the external potential. Finally, all terms of the Euler equation Eq. (pS)) 
are lead back to a linear density fluctuation 5nc{r^t). Thus, one ends up with 



4 



ern ^ — ern 



2A / V 2A 



[l+e^"=«"]erf (':|) ). (32) 



This relation leads to real valued solutions for the resonance frequencies for standing waves with fc„ — wk / only 
and scales with the plasma frequency Wpi. In the limit fc — > we find w(0) = Wpi which coincides with the bulk limit. 

From the eigenvector of the plane wave mode, one can derive the wave number k — 'k / Rc, which corresponds to 
n — 1. This means the dispersion of the plane wave mode is determined by the radius of the electron cloud. Results 
for this case are shown in Fig. [S] (b) for different cluster sizes and are compared with the simulation data. For the 
cluster with 1000 ions a plane wave mode with k = 37r/i?e and n = 3 was found as well. In Fig. |6] (b), it is marked 
with an empty square. Its spatial structure is shown in Fig. [S] (d). The simulation data for the Naiooo cluster fit the 
dispersion Eq. p2[) as well. Deviations of the plane wave resonance for smaller clusters are caused by to the radial 
dependence of the electron density profile. 
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VI. CONCLUSION 

We have investigated collective excitation modes of a nano plasma in highly excited metal clusters. The collective 
excitation of electrons inside the cluster are obtained from bi-local current density correlation functions by solving 
the eigenvalue problem of the current-density correlation matrix. Using RMD simulations, the local current density 
j{f, t) for excited clusters of 55 up to 1000 ions with densities of n[ = 2.15 • 10^^ cm~'^ as well as 2.80 • 10^^ cm~^ and 
temperatures of Tc = 1 eV have been investigated. Pseudo-potentials of sodium were used to calculate the electron 
dynamics without consideration of degeneration effects any further. For the analysis of electron dynamics at lower 
temperatures, the inclusion of quantum effects for the calculation of the local current density j{r, t) of cluster electrons 
is an open question at this point. It would be useful to go beyond present classical description to discuss for example 
cold, non-excited clusters. 

The spectrum of dipole-like modes was investigated in more detail. Using analytical calculations, it was possible 
to relate the position of resonance modes in the frequency domain to their spatial mode structure. Results for the 
cluster size dependence of the resonance frequency have been shown. A smooth transition to the bulk behavior has 
been obtained. The analysis of further resonance frequencies and also other modes including breathing modes would 
be desirable. The width of mode resonances and the role of collision-less damping effects as well as the collision 
frequency need to be investigated in the future. The systematic change of the collision frequency with cluster size up 
to the bulk limit remains an interesting field. 

From RMD simulations, different collective excitations have been found in nano plasmas, including dipole-like 
and breathing modes. These collective excitations will influence the scattering and absorption properties of clusters, 
see [l3|. Collective effects of electron motion play a role when analyzing ultraviolet (UPS) or x-ray photo-electron 
spectroscopy (XPS) experiments, as has been pointed out by Andersson et al. [22|. It is a challenge to experimentalists 
to confirm the occurence of different collective excitations in nano plasmas. 
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